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Abstract. - We numerically study the interplay between heterogeneous dynamics and prop- 
erties of negatively curved regions of the potential energy surface in a model glassy system. 
We find that the unstable modes of saddles and quasi-saddles undergo a localization transition 
close to the Mode-Coupling critical temperature. We also find evidence of a positive spatial 
correlation between clusters of particles having large displacements in the unstable modes and 
dynamical heterogeneities. 



The dynamics of supercooled liquids is often described as a complex trajectory across their 
Potential Energy Surface (PES) [1]. This approach traces its origin back to the pioneering 
work of Goldstein [2], who argued that at low temperature these systems get trapped in the 
basins of attraction of local minima of the PES and hence the atomic dynamics is slowed 
down. In the last decade, several authors have addressed a quantitative study of the PES 
through numerical simulations of model systems. The original picture has evolved into a 
refined description of the dynamics in terms of collections of local minima (metabasins [3]) 
and transitions between them [4-6] . 

Other works have focused on the properties of the Hessian matrix of the potential energy, 
revealing even a more complex scenario, where not only local minima but also higher order 
stationary points (saddles) and more general points characterizing regions of negative cur- 
vature of PES (quasi-saddles) play an important role in the structural slowing down of the 
liquid [4-15]. Indeed, it has been shown [12,13] that some information about the liquid-like 
diffusive dynamics is encoded in the imaginary spectrum of the Hessian matrix: the num- 
ber of unstable modes n im of saddles is correlated with the diffusivity of supercooled model 
systems and decreases as the liquid is cooled. Attempts have been made [10-13] to identify 
the temperature at which the thermal average of n im extrapolates to zero with the critical 
temperature T c where the purely dynamical Mode-Coupling Theory (MCT) [16, 17] predicts 
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a structural arrest of the liquid. The residual relaxation exhibited by both real and simulated 
systems below T c is to be attributed, within this scenario, to "activated processes" , i.e. rare 
transitions over finite energy barriers, ignored in the mean field approach of MCT. 

However, the idea that the slowing down of the liquid around T c reflects a geometric transi- 
tion [11] in the PES has been recently criticized. Careful studies have in fact put in doubt the 
existence of a sudden change in the sampling of saddles at T c [4, 15] and there are indications 
that saddles do not disappear completely below T c [14,18]. Most worrying, perhaps, is that a 
saddle-based approach appears unsuitable [19] for the description of dynamical heterogeneities 
[20-22], i.e. rearrangements involving localized subsets of mobile particles, which have been 
recognized in the last years as the hallmark of the supercooled dynamics and whose evidence 
is hard to glimpse on a 3N-dimensional PES. 

In this letter, we analyze the localization properties of the unstable modes of saddles and 
quasi-saddles sampled by a supercooled Lennard-Jones (LJ) mixture. We provide quantitative 
evidence that the spatial structure of these unstable modes change from extended to localized 
around T c . We establish for the first time the existence of a direct connection between saddles 
and dynamical heterogeneities in a supercooled liquid. In fact, the regions of the system where 
unstable modes are localized are statistically correlated with clusters of particles which are 
highly mobile on the /3-relaxation timescale. Our results thus indicate that a crossover in the 
nature of the sampled saddles is indeed underlying the emergence of dynamical heterogeneities 
and the failure of mean field theories below T c . 

Our model is the binary LJ mixture introduced by Wahnstrom in [23]. This mixture has 
been shown to be a good glass-former and its properties have been analyzed in-depth [9,23], 
especially as far as dynamical heterogeneities are concerned [24] . Our system is an equimolar 
mixture of 500 particles interacting via the LJ potential u a fj{r) — 4e Qi g[(cr Q ^/r) 12 — (cr Qi g/r) 6 ], 
with a, (3 = 1,2 indexes of species, and enclosed in a cubic box of side L = 7.2779 an with 
periodic boundary conditions. Reduced units will be used in the following, assuming an, 
en, and (mur^/en) 1 / 2 respectively as unit of distance, energy and time. The interaction 
parameters are on = 1.0, a 2 2 = 0.833, aw = 0.917, en = e 2 2 = e i2 = 1.0. The masses are 
mi = 1.0, m2 = 0.5. Two different cut-off schemes for the potential have been considered, 
cut and shifted (CS) and cut and quadratically shifted (QS) at r c = 3.0. The QS cut-off 
is obtained adding a term a + br 2 with a and b determined to ensure continuity to both 
u(r) and its first derivatives at r c . This cut-off scheme has been used to check any bias in 
the determination of the properties of the PES due to discontinuities in the force at the cut- 
off [25]. Molecular Dynamics (MD) simulations have been performed in the canonical ensemble 
using the Nose-Poincare thermostat [26,27]. At the lowest temperature corresponding to an 
equilibrated supercooled liquid (T = 0.575), the length of the run was 2 x 10 7 time steps 
with St = 0.008, equivalent to a total simulation time two orders of magnitude larger than 
the structural relaxation time. Our best estimate of the MCT critical temperature, obtained 
from an analysis of the intermediate scattering functions, is T c = 0.55 ± 0.01, which is slightly 
lower than the previously reported ones [9]. 

From each equilibrated run we extracted between 100 and 200 independent configurations 
and minimized the mean squared force W = ^/NJ2i ff using the LBFGS algorithm [28] to 
locate the nearest saddle. It has been recognized [12,14,25] that ^-minimizations usually stop 
near quasi-saddles, i.e. local minima of W where the PES displays one inflection mode with 
non-zero force, and that discontinuities in the cut-off can lead to incomplete minimizations [25]. 
While other numerical strategies have been proposed for accurately locating saddles [15], there 
has also been evidence [12,13,29] that quasi-saddles do not differ significantly from true saddles 
as far as the influence on the dynamics is concerned. In the following we will use the term 
saddle in a wide sense, without distinction between quasi-saddles and saddles. The Hessian 
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Fig. 1 - Gyration radius L^*" (circles, left axes) and participation ratio p %m (squares, right axes) as 
a function of T, both in the supercooled liquid (black) and glass (white). Error bars represent one 
standard deviation of the distribution of values. The horizontal dashed line is drawn at Lt ox /2. 



matrix at each saddle was diagonalized yielding the eigenvalues v 2 , and the eigenvectors e", 
where a = 1, . . . , 3N is an index of mode and i = 1, . . . , N an index of particle. The order 
rii m of a saddle is then defined as the number of imaginary frequencies in its spectrum. 

To have a better understanding of the localization properties of the unstable modes, we 
consider the average 



1 



a = l 



(1) 



of the squared displacements ef 2 of atom i over the rii m unstable modes. The vector E" n = 
(El" 1 , . . . , E 1 ^ 1 ) contains an averaged information about the distribution in real space of the 
instabilities associated with the saddle. By construction E" n is normalized so that we can 
define a participation ratio for the unstable modes of a given saddle in the usual way [26] 



= N 




(2) 



p ml will be roughly 1 for extended instabilities and 0(1/ 'N) for a single, localized one. Further 
insight can be gained considering the gyration radius [30] 
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where r g = J2i fiE™ 1 . For extended instabilities L l c 



c /2. Both p im and Lf 1 refer to a 



given configuration and quantify the degree of localization of the unstable modes of a saddle 
as a whole. They will overestimate the size of localized instabilities when the system is large 
enough that several independent ones are present, but we will show that this should not be 
the case for the mixture in consideration. 

We show in fig. ^the thermal average of the participation ratio p %m and gyration radius 
-L™ as a function of T. At high T saddles are extended: p lm tends to saturate around 0.7 and 
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Fig. 2 - Plot a) : distribution of average displacements E\ m on the unstable modes of a saddle sampled 
at T = 0.575. Plot b): normalized propensity of motion Ef of the MD configuration from which the 
saddle in plot a) has been located. A sphere proportional to E\ m and E\ is drawn around each 
particle, in plot a) and b) respectively. 

L 1 ™ approaches its limiting value of Lb ox /2. In the landscape-influenced range [8], roughly 
between T c and 2T C , the instabilities of saddles have a mixed character. In fact, while p lm 
decrease markedly, the gyration radius L*™ still fluctuates around L{, ox /2, indicating that 
most of the instabilities of saddles are percolating through the system despite their enhanced 
localization. On approaching T c saddles become strongly localized and the average L* m de- 
creases abruptly. Figure^thus show that the localized nature of saddles becomes pronounced 
exactly in the same temperature range where dynamical heterogeneities first appear [19] and 
in correspondence to the break-down of mean-field approaches. 

The degree of localization of the instabilities of saddles can be inspected directly in the 
snapshot of fig. which displays the spatial distribution of E\ m in the simulation box for a 
saddle sampled at T = 0.575, slightly above our estimated T c . A sphere of radius proportional 
to E\ m is drawn around each particle. We see that particles involved in the unstable modes, 
i.e. those with a large E\ m , are strongly clustered. As counterparts, extended regions of 
stable particles are observed. The spatial structure of saddles like that in fig. [2^, for which 
n-im = 7, casts some new light on a common assumption about the nature of saddles sampled 
by supercooled liquids, namely that saddles of low order may originate from non- interacting 
subsystems, each experiencing a saddle of order one [13, 15,31]. Figure^ shows that the 7 
unstable modes have a strong spatial overlap, being concentrated on a rather compact cluster 
of particles. Similar features are present in the whole landscape-influenced range, while in the 
hot liquid the fraction of extended unstable modes increases so that particles are essentially 
involved in all of them. Thus, the independent subsystems interpretation does not seem to 
hold for saddles sampled by LJ systems in the T > T c range ( 1 ). 

The emerging scenario seems to mirror the observations of dynamical heterogeneities [19]. 
The question arises naturally whether there is a direct mapping between saddles and dynamical 

( 1 )Similar results were in fact obtained for the BMLJ of [35]. We have also checked that this is the case also 
for the true saddles we have sampled, for which W f» 10~ 10 . 
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heterogeneities, i.e. whether unstable clusters like that in fig. are also the mobile ones. 

To address this issue we followed the approach of [32], who introduced the propensity of 
motion of particles as a measure of dynamical heterogeneities. A configuration is selected 
from the MD trajectory at a given temperature T. Then several short runs at the same T 
are performed, monitoring the square displacements Arf(t) from the reference configuration. 
We evaluate the propensity of motion as (Ar?(i*)), where (...) indicates an average over 
independent sets of initial velocities and t* corresponds to the maximum of the non-Gaussian 
parameter a-i (t) , which lies in the late /3-rclaxation regime and is a characteristic timescale for 
dynamical heterogeneities [20]. Using this procedure it is possible to identify the dynamical 
heterogeneities that a single configuration will, on average, give rise to. We have considered 
up to 500 sets of initial velocities, drawn from the appropriate Maxwellian distribution. Since 
the timescales for diffusion are slightly different for the two species (roughly a factor of 2), 
we have also tried to choose different values of t according to the species, but the picture was 
essentially unaltered, so we simply used the value of t* for the small particles. 

Figure [5J} shows the distribution in the simulation box of the normalized propensity of 
motion 

«E£i(Ar?(**)) a » 1/2 

for the MD configuration whose nearby saddle is actually that of fig. [3a. By comparison 
we see that the localization and essential morphology of the mobile cluster identified by the 
propensity of motion of the MD configuration are well reproduced in the cluster of average 
unstable displacements of the nearby saddle. To assess the statistical relevance of the cor- 
relation we have analyzed 20 independent configurations at T— 0.575. We define as mobile 
(immobile) those particles for which the normalized propensity of motion is larger (smaller) 
than a threshold e h (e l ). Applying the same cut-off procedure to the average displacements 
in the unstable modes E\ m , we introduce an analogous separation in unstable and stable par- 
ticles. The two threshold values e h and e l could vary or even coincide without altering the 
overall picture. With our choice e l = 0.01, e h = 0.045, the fraction of particles belonging on 
average to each subpopulation is around 25%. 

In fig. |2| we show the radial distribution functions at T— 0.575 for mobile-unstable (MU), 
immobile- stable (IS) and mobile-stable (MS) pairs, compared with the total radial distribu- 
tion function g(r) = (gn(r) + 2t/ 12 (r) + <?22( r ))/4- The significant enhancement of the first 
two coordination shells in gMu( r ) an d gis( r ) clearly shows that particles with a high (low) 
propensity of motion are surrounded, on average, by particles having large (small) displace- 
ments in the unstable modes of the closest saddle. Consistently with such correlation, the 
9Ms{ r ) (identical to gsAi(r)) is lower than the total g(r). Moreover, as expected, the effect 
decreases by increasing temperature. All this provides statistical evidence that dynamical 
heterogeneities have spatial correlation with compact clusters of particles taking part in the 
unstable modes of nearby saddles. Interestingly, fig. also suggests that regions of immobile 
particles possess distinct structural properties, as evidenced by the deeper separation between 
the first and second shell of neighbours and the enhanced splitting of the second peak in gis(r) 
and in agreement with [21]. 

In summary, by performing MD simulations for a supercooled LJ mixture we have shown 
that the relaxation channels associated with the unstable modes of saddles and quasi-saddles 
sampled along the MD trajectory crossover from spatially extended to localized around T c . 
These localized unstable modes display non-trivial spatial correlations with the dynamical 
heterogeneities identified by the propensity of motion [32] in the /3-relaxation timescale. The 
novel finding of such a correlation may represent a step ahead in the understanding of the 



6 



EUROPHYSICS LETTERS 



12 




r 



Fig. 3 - From top to bottom: radial distribution functions for mobile- unstable (MU), immobile-stable 
(IS) and mobile-stable (MS) pairs of particles, compared with the total distribution function (dashed 
lines) at T = 0.575 (see text for definitions). Data related to the IS and MU cases have been shifted 
for clarity. Error bars represent one standard deviation. 

heterogeneous nature of the supercooled dynamics. Theoretical and computational tools that 
have been developed for the study of the PES [33,34], may now reveal their utility in a direct 
analysis of dynamical heterogeneities. We note that other observablcs, e.g. free volume, have 
been recently addressed as possible origins of dynamical heterogeneities, but their spatial 
correlation with local mobility has been cither found to be poor [36] or investigated below 
the glass transition [37]. In this perspective, our results point to a key role of saddles and 
quasi-saddles in bridging the gap between the PES description of supercooled liquids and 
approaches which focus more directly on the dynamics like the heterogeneity-based picture 
proposed by Bcrthier et al. [19] or the mean field MCT. 

* * * 

We would like to thank F. Sciortino for useful comments. 
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